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Abstract 

We compute the neutralino relic density in the minimal supersymmetric standard 
model by using exact expressions for the neutralino annihilation cross section into 
all tree-level final states, including all contributions and interference terms. We find 
that several final states may give comparable contributions to the relic density, which 
illustrates the importance of performing a complete calculation. We compare the 
exact results with those of the usual expansion method and demonstrate a sizeable 
discrepancy (of more than 10%) over a significant range of the neutralino mass of up 
to several tens of GeV which is caused by the presence of resonances and new final- 
state thresholds. We perform several related checks and comparisons. In particular, 
we find that the often employed approximate iterative procedure of computing the 
neutralino freeze-out temperature gives generally very accurate results, except when 
the expansion method is used near resonances and thresholds. 
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1 Introduction 



One of the most exciting unresolved issues in particle physics and cosmology is the 
nature of dark matter (DM) in the Universe. Cosmological and astronomical obser- 
vations have put, over the years, significant restrictions on its expected properties. In 
particular, it is widely expected that a hypothetical candidate should be in the form 
of cold DM (CDM), while a combination of several measurements has narrowed the 
range of its relic abundance to 

0.1 <n CD Mh 2 £0.3, (1) 

where ^cdm is the ratio of the CDM relic density to the critical density and h ~ 0.7 is 
a parameter in the Hubble constant H = 100 h km/ sec/Mpc Q. Constraints from the 
Big Bang Nucleosynthesis || and, recently, from observations of the cosmic microwave 
background radiation j3],|J have strongly restricted the abundance of baryonic DM; 
for example, Ref. || quotes the range Qf,h 2 = 0.02 ±0.002. Therefore most of the DM 
in the Universe is expected to be non-baryonic. 

From a particle physics point of view, CDM is most naturally made of some weakly 
interacting (stable) massive particles (WIMPs). Supersymmetry (SUSY) ||, which 
over the years has emerged as a leading candidate for new physics beyond the Standard 
Model (SM), predicts that, in the presence of R-parity, the lightest supersymmetric 
particle (LSP) is stable. This alone makes the LSP an interesting candidate for a 
WIMP and DM in the Universe. 

A SUSY candidate for the LSP WIMP is by no means unique. In models coupled 
to gravity there exists the gravitino, the SUSY partner of the graviton. The gravitino 
has long been known to be a potential candidate for DM, although it generically 
suffers from a well-known "gravitino problem" |7]]. In SUSY models that incorporate 
the Peccei-Quinn solution to the strong CP problem, there exists the axino, the 
fermionic partner of the axion. The axino has recently been shown to be an attractive 
and well-motivated alternative candidate for the LSP and CDM ||. A number of 
more speculative possibilities, stemming from, for example, string theories, have also 
been proposed. However, in some sense the simplest choice is to consider one of 
the superpartners of the SM particles as a potential LSP. Among these, the lightest 
neutralino \ stands out as probably the most natural and attractive candidate for the 
LSP and DM @,|10[- In this paper we will focus on this case. 
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In order to be able to perform a reliable comparison between theoretical predic- 
tions and improving measurements of the relic abundance and limits from under- 
ground DM searches, a precise calculation of the relic density becomes necessary. The 
constraint ([!]) is known to often provide a strong restriction on allowed combinations 
of SUSY masses and couplings. For example, in the minimal supergravity model (also 
dubbed CMSSM) the condition fl x h 2 < 0.3 provides a strong upper bound of a few 
hundred GeV on the masses of gauginos and scalars over an overwhelming fraction 



of the parameter space [p],|Tz| . When combined with recent LEP and other data, the 
remaining CMSSM parameter space becomes highly restricted ||13|| . In the future, the 
range (|1|) is expected to be significantly narrowed by MAP and Planck, which in turn 
will lead to even stronger constraints on allowed combinations of SUSY masses and 
couplings. 

One of the key steps in computing the relic density of the neutralino involves 
calculating the thermal average of its annihilation cross section times the relative 
velocity. This calculation is technically rather involved and, therefore, in most of 
the literature, approximate methods have been used. These methods are based on 
expanding the thermal average in powers of x = T/m x , where T is the temperature of 
the thermal bath and m x is the neutralino mass. One usually computes only the first 
two terms of the expansion. It has been well known that the expansion fails badly near 



s-channel resonances [ 14 — IT | and near new final-state thresholds [14, 15 1 . In particular, 



it was emphasized in Ref. |L6|] that, owing to the very narrow width of the lightest 
SUSY Higgs boson h, the error can be as large as a few orders of magnitude very close 
to the /z-resonance. However, even though the usual expansion is generally expected 
to be accurate enough further away from resonances and thresholds, no detailed study 
of this point has yet been done. 

On the other hand, while a general formalism for computing the thermal average 
precisely was derived a long time ago [T^,[TB[, a full set of exact, general and explicit 
analytic expressions for the neutralino pair-annihilation cross section is still not avail- 
able in the literature. (The thermal average can be obtained by performing a single 
integral over the cross section, as we will see below.) The cross sections for the neu- 
tralino pair-annihilation into the SM fermion-pair (//) and the lightest Higgs-boson 



pair (hh) final states were given in Ref. fl7 [, and for the WW and ZZ final states 
in Ref. ||19|| . The annihilation into // is often most important. However, other final 
states can also play a considerable role, depending on the model. 
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We have derived a full set of exact analytic expressions for the cross section of the 
neutralino pair-annihilation in the minimal supersymmetric standard model (MSSM). 
We have made no simplifying assumptions about the degeneracy of the left- and right- 
sfermion masses, and we included all tree-level final states and all intermediate states. 
We have kept finite widths in s-channel resonances. While interference terms and 
usually sub-dominant channels are sometimes neglected in the literature, in some 
cases they may give significant contributions. In particular, the channels involving 
gauge-Higgs boson may in some cases be sizeable or even dominant. One example 
where this is the case is the W ± H Zf final state, and in this paper we give expressions 
for the cross section for the chargino exchange contribution. A full set of all the 
expressions for all the final states will be presented elsewhere 

As a check, we have performed a numerical comparison of our cross section with 
the results obtained by using the recently released code DarkSusy pl| . We have 
found, for the same values of input parameters, an impressive agreement, at the level 
of a few per cent, for all the annihilation channels, which we find reassuring. 

We have also made a comparison of the relic abundance computed using our exact 
formulae with the one obtained using our expansion formulae. So far, no such detailed 



analysis has been presented in the literature in the case of the MSSM. In Ref. [17 



an analogous comparison was made in the context of the highly constrained minimal 
supergravity scheme and only relatively close to resonances. We confirm that the 
expansion gives highly inaccurate results near resonances and new thresholds. We 
also show that very far from such cases the error is typically rather small, of the 
order of a few per cent. However, we find that, because of the existence of several 
resonances (Z and the Higgs bosons), the expansion produces large errors, compared 
to an exact treatment, over a sizeable range of m x , even of a several tens of GeV. We 
therefore conclude that the widely used method of expansion may lead to significant 
errors in a sizeable fraction of the neutralino mass which is typically expected not to 
significantly exceed ~ 200 GeV by various naturalness criteria . 



In the current version of our code we have used a popular approximate iterative 
method of determining the freeze-out point. We have examined the accuracy of the 
procedure and compared our results with the ones obtained using DarkSusy. We have 
found that the method in general works very well when both exact and expansion 
methods are used, except in the latter case near resonances and thresholds where it 
may even break down. The iterative procedure can be safely applied in such special 
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cases, provided exact expressions for the cross section are used and much care is paid 
to properly integrating the thermal average numerically. 



2 Calculation of the Relic Density 

The relic abundance of some stable species x is given by ft x = p x j p cr it, where p x is 
the relic's (mass) density and p cr u = 3Hq/8tiG = 1.9 x 1CT 29 (h 2 ) g/cm? is the critical 
density. It can be computed by solving the Boltzmann equation, which describes the 
time evolution of the co-moving number density n x in the expanding Universe,[] 

^ = -3Hn x - (av M0l ) (n\ - (n?) 2 ) , (2) 

where n x q is the number density that the species would have in thermal equilibrium, 
H(T) is the Hubble expansion rate, <j(xx ~ > a U) denotes the cross section of the 



species annihilation into ordinary particles, vm_0\ is a so-called M0ller velocity |15 
which is the relative velocity of the annihilating particles, and (crvu0i) represents the 
thermal average of uvm0\ which will be given below. In the early Universe, the species 
X were initially in thermal equilibrium, n x = n x q . When their typical interaction rate 
T x became less than the Hubble parameter, T x < H, the annihilation process froze 
out. Since then their number in a co- moving volume has remained basically constant. 

In order to calculate the relic density precisely enough, one has to consider the fol- 
lowing ingredients: (i) compute all the contributions to the annihilation cross section, 
(ii) use an exact formula for the thermal average, (iii) rigorously solve the Boltzmann 
equation, and (iv) include co-annihilation. (This becomes important when the mass 



of the next LSP is nearly degenerate with the LSP mass |14| , |26| - [2ql .) In this work, we 
will concentrate on points (i) and (ii) but will also examine to some extent point (iii). 
In particular, we will compare the usually used expansion formulae with the exact 
ones. We will not consider here the effect of co-annihilation. 

An often used, approximate, although in general quite accurate (see later), solution 
to the Boltzmann equation is based on solving iteratively the equation 



x 



' - MSJ^<<"'m»>>(*/)4 /2 ). (3) 



2vr 3 V 2g,G 



N 



1 For a review of calculations of the relic density, see, e.g., Refs. 
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where represents the effective number of degrees of freedom at freeze-out (y^T — 9) 
and Gn is the Newton constant. Typically one finds that the freeze-out point Xf = 
Tf/m x is roughly given by 1/20. The relic density p x = m x n x at present is given by 

1.66 /T y \ 3 „ ^_ 1 

where Mpi denotes the Planck mass, T x and T 7 are the present temperatures of the 
neutralino and the photon, respectively. The suppression factor (T x /T 7 ) 3 ps 1/20 
follows from entropy conservation in a comoving volume | 29fl . Finally, J{xj) is given 
by 



J(xt) = / dx(av M 0i)(x), (5) 

JO 

where x = T/m x , as defined earlier. 

Below we will concentrate on computing J(xf) and in particular on the thermally 
averaged annihilation cross section times velocity (<7t>M0i)- A proper definition which 
includes separate thermal baths for both annihilating particles is given by |]TS|,0 



, ™ _ fd 3 Pl d 3 p 2 av M0l e- E ^ T e- E ^ T 

{aVM0l} (T) " Jd^p 2 e-^e-^ (6) 

where p\ = (E\, pi) andp2 = (E 2 , p 2 ) are the 4-momenta of the two colliding particles, 
and T is the temperature of the bath.Q The above expression can be reduced to a 
one-dimensional integral which can be written in a Lorentz-invariant form as |15| 



(av M0l )(T) = 1 f°° dsa(s)(s - 4m x )^i (ffi , (7) 

8m x TK 2 (m x /T) J 4m| x \ T J 

where s = (pi + p 2 ) 2 is a usual Mandelstam variable and Ki denotes the modified 
Bessel function of order i. 

In computing the cross section <j(s), it is convenient to introduce a Lorentz- 
invariant function w(s) [0 

w(s) = ^JdUPS\A(xX^&\\)\ 2 (8) 

where |.4.(xx — > all) | 2 denotes the absolute square of the reduced matrix element for 
the annihilation of two x particles, averaged over initial spins and summed over final 



2 Note that in many cases one uses another definition of < avuai > which involves a single thermal 
bath for both neutralinos. Compare, e.g., Refs. p4,E3,p4|. 
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spins. The function w(s) is related to the annihilation cross section via [17 



For a general annihilation process into a two-body final state XX 
function w(s) can be written as 



wis) 



32 7i 



Yl e { s ~ ( m /i + m / 2 ) 2 ) Pf(s,m fl ,m f2 )w hf2 (s) 



where the summation goes over all possible channels, and 



W hh{ S ) 



-J- fdn\A(xx^hf2 

07T J 



where 



P f (s,m fl ,m h 



{m fl +m h y 



1/2 r 



1/2 



(9) 

/1/2, the 
(10) 

(11) 
(12) 



Computation of w(s) is in general rather involved and final analytic expressions 
exhibit considerable complexity, as will be illustrated in the next Section. This is one 
reason why in most of the literature one uses the well-known approximation in terms 
of expansion in powers of x, (ctvm0\) — a + bx. Using the definition (|6|), the expansion 



reads 18 



< crv M 0\ >= — r 



w - - (2w-w')x + 0(x 2 ) 



a + bx + 0{x 2 



(13) 



where w'(s) denotes dw{s)/d{s/Am^}. 

Note that in case of the expansion of < &VM01 > defined in terms of a single thermal 
bath for both neutralinos, the corresponding coefficient a' is the same as Eq. flllf ) while 
b' = b+ fa g. For more details, see, e.g., Ref. plj . 

The expansion method is widely used not only because it is computationally some- 
what less involved but also because it is expected to be relatively accurate. This can 
seen by examining the behavior of the integrand in Eq. ([?]). For a massive particle like 
the neutralino for which T < m x /20 and ^/s > 2m x , the argument of the function K\ 
is much larger than unity. Since K\{y) ~ ^Ti/2ye~ v at y ^> 1, the thermal average (0) 
can be written as a convolution of the cross section with a function f(s): 



(o-v M0 i) 



4m| 



dsa(s)f(s) 



(14) 
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where f(s) has an exponential suppression factor, f(s) oc e v^/ T _ Thus one would 
expect that the usual expansion in powers of x should converge quickly [[i~7j . 



One should however use this argument with some care. First, the annihilation 
cross section may change rapidly with s. It is well-known that this happens, e.g., 
near s-channel resonances and thresholds of new final states. In such cases the expan- 



sion method produces large errors |T4]-|T7||. Far away from such singular points the 



expansion is expected to give accurate results but, to our knowledge, this has never 
been explicitly verified in the literature. Finally, one would want to examine actually 
how far from resonances one can start relying on the usual approximation. As we will 
see later, the range of m x where the expansion fails to be reasonably accurate can be 
actually quite large. 

3 Analytic Results 

As stated in the Introduction, we have derived a full set of exact, analytic expressions 
for the cross sections for the neutralino pair-annihilation processes into all allowed 
(tree-level) final states in the general MSSM. We have included all contributing di- 
agrams as well as all interference terms and kept finite widths of all s-channel reso- 
nances. We have made no simplifying assumptions about sfermion masses although 
we assumed that there are no CP violating phases in SUSY parameters. Here we will 
give one analytic example of our exact expressions. A full set of results for all the 



final states will be given elsewhere f20| 



Next, starting from the exact expressions we have computed the coefficients a 
and b using Eq. (|T3|) for all the partial annihilation channels. In the literature there 



exist several analytic formulae for the expansion coefficients, including [p5|,|32|-p4 
but, due to different conventions and complexity of expressions, comparison is not 
always doable. We have checked our results for the a-coefncients in appropriate limits 
against published results and agreed in all cases. In the next Section we will present 
a numerical comparison of the approximation in terms of the usual expansion with 
our exact results. 

Let us begin by introducing the relevant MSSM parameters.^ The lightest neu- 
tralino is a mass eigenstate given by a linear combination of two neutral gauginos and 



3 We follow the convention of Ref. 
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two neutral higgsinos 

X = N n B + N 12 W 3 + N 13 H° + N 14 Hl (15) 

The neutralino mass matrix is determined by the U(l)y and SU(2)y gaugino mass 
parameters Mi and M 2 (and we impose the usual GUT relation M l = § tan 2 6 W M : 
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the Higgs/higgsino mass parameter \x and tan/5 = v 2 /v\ which is the ratio of the 
vacuum expectation values of the two neutral Higgs fields. 

There are two neutral scalar Higgs bosons h and H and one pseudoscalar A plus 
a pair of charged Higgs H . (We will typically suppress the Higgs charge assignment 



except where this may lead to ambiguities.) We use expressions given in Ref. |36| for 
computing radiatively corrected Higgs masses. 

Other relevant parameters which determine the masses of scalars and various cou- 
plings are the squark soft mass parameters rriQ, my and mp, the slepton soft mass 
parameters my and m#, and the pseudoscalar mass m^. We also include the trilinear 
terms A% (i = t, b, r) of the third generation which are important in determining the 
masses and couplings for the stop, sbottom and stau sfermions, respectively. 

There are a number of final states into which the neutralino can pair-annihilate. 
They include: X X -> //, WW, ZZ, Zh, ZH, ZA, W ± H^, AA, AH, Ah, HH, Hh, 
hh and H + H~ . Among these, in the gaugino region, fermion pair // final states 
usually give dominant contributions, unless the exchanged sfermion masses mj are 
large and in some special discussed earlier. These final states are also always 

kinematically allowed (except for tt) but, due to the s-wave suppression 0, their cross 
section is proportional to m^/m^r. For the higgsino-like neutralinos the states ZZ and 
WW become very important once kinematically allowed ||19| . However, these final 
states can give significant contributions even in the case of the gaugino-like neutralino. 
This is because WW and ZZ are not s-wave suppressed, unlike the channels hh, Hh, 
HH, AA, H + H~ and ZA | ]25[| . What is interesting is that the cross sections for the 
other s-wave unsuppressed states W ± H T , Zh, ZH, Ah and AH, once kinematically 
allowed, can be even larger than those of WW and ZZ. In fact, they can even 
dominate over those of //. These points will be illustrated with numerical examples 
in Section [|[ 

Here, we present a set of expressions for a t- and ^-channel chargino (X12) exchange 
contribution to the process XX ~ > W + H~. This contribution is often dominant. 
(There are also s-channel diagrams involving h, H and A exchange. Since typically 
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ruh "C ran ~ ~ ttih± } the Higgs resonances will be outside of the kinematically 
allowed region.) The function w, Eq. ([II]) in this case reads 



w 



w+h- 



rn 



E 

W i,j=l 



r(2) , rC 1 ) i r(°) 



(16) 



where 

,-(2) 

ij 



7^ - (s - - m^± - 2 m£) 71 + % + 6 m^m^ J 



+:i 



6 m^m^ %-y 2 + g^ y 



(17) 



r(l) 



Re(C5 w *C?f + D 



HW* t^HW^ 

+i u -j , 



AO) 

ij 



2T 2 + 2(2m 2 x -m%y)Ti 
- 2{m 2 x - m 2 w ){m 2 x + 2 m 2 w )% + 2 m 2 w y x - y ' 

+ B*(C*TC™ - DirV™) Kr 1 + 6 m 2 w (m 2 x - m 2 H± )T 

-2y 2 -( s -m 2 H± -4 m^ r )y 1 - y } , (18) 

{C hw* c hw + d^TdBW) [( s - m 2 x -2 m 2 w ) T 2 - G%% T x 

- G^% + {s-2m 2 v )y 2 + G^y ' 



HW* nW^ 



D i D 



6 m 2 w m 2 Ti - m 2 y 2 + y 



(19) 



Obviously the total contribution to w(s) from W ± H Zf is twice that from W + H~. 
The symbols and in the above equations are combinations of several 

coupling constants: 



r*HW _ n x i xH n x t xW _i_ r> x T xH n x t xW ~ 

L/ ±i = °S °F =■= °P U A 



(20) 
(21) 



The relevant couplings describe neutralino-chargino-VT / H interaction terms in the 
Lagrangian as 
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+ E XT {cf XH ~ - C£ XH ~ 15 ) X H- h h.c. (22) 



in the convention of Ref. |35|. The two pairs of charginos are denoted by xf (i — 1 ; 2) 



and m x ± represent their masses. The functions T k = T k (s,m^,m^ ± ,m^,m^,m^±), 

34 = yk{s,ml,mH±,m^,mK,mK), k = 0,1,2, and the functions , G$]f 5 > 

are listed in the Appendix. 

4 Numerical Results 

In this Section we will illustrate the points made above by presenting several numerical 
results. We will concentrate on the gaugino-like neutralino. We remind the reader 
that a nearly pure gaugino (bino) is the most natural choice for the CDM both in the 
MSSM |22| and in the CMSSM where in most cases the LSP comes out to be a nearly 
pure bino ||iT|,|T^,[T6| . 



We begin by plotting in Figs. [I] and [| the function J{xf), Eq. (H), computed 
using exact and expansion methods as a function of m x . We take tan/3 = 10, 
m = rag = ra v = m D = m L = m E = 500 GeV, = 400 GeV, A t = A b = 600 GeV 
and \i = —500 GeV. To obtain sparticle and Higgs mass spectra we have used the 
code DarkSusy. While several analytic expressions and numerical codes (e.g., SUS- 
PECT |J7|]) are available in the literature, here we have used DarkSusy to facilitate the 
comparison of our results for the cross sections with those produced by the package. 
For the values given above we obtain the following Higgs masses: = 117.3 GeV, 
m H = 400.5 GeV and m H ± = 407.9 GeV. 

We show both the total and the individual contributions from all the final states. 
The solid and dotted lines represent the exact and the expansion results, respectively. 
For the sake of comparison of both methods, for now we do not impose experimental 
constraints. Note that m x increases with increasing M2 while we keep \i fixed at a 
rather large value. This means that in these Figures (and also the following ones 
below), the neutralino is mostly a gaugino. Even at m x = 400 GeV, + Nf 2 > 0.9. 
(However, between 400 GeV and 500 GeV the LSP turns into a nearly pure higgsino.) 
In this region the effect of co-annihilation with the lightest chargino and the next 
lightest neutralino, which we do not consider here, is not important and can be safely 
neglected. 
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For each final state, we compute J(xf), Eq. @, numerically via Eq. (|7|) by using 
our exact expressions for the cross sections ( Jexact (£/))• Next, we use the cross sections 
to compute the first two expansion coefficients a and b, as in Eq. (|13D, and to compute 
Jexp(xf) = axf + \bx 2 f. For consistency, in both cases we determine the freeze-out 
point xj by employing the iterative procedure (|3|). We will come back to the accuracy 
of the iterative procedure below. 

In the behavior of the total J{xf) one can clearly recognize three peaks which are 
due to the s-channel Z, h and A/H exchange contributions. In particular, we can 
see that the A-resonance is very broad. (Due to a near mass degeneracy itlh — niA, 
the usually narrower if-pole is 'burried" underneath the v4-pole.) This is caused by 
the amplification of the coupling Abb ~ tan/3 at larger tan/? and the fact that, in 
the resonance region, the amplitude- square for the A-exchange is proportional to s, 
unlike the case of h and H where it goes like s — 4m^ « at ^/s ~ 2m x . The Z-pole 
is also clearly visible in all the partial // final states while the Higgs poles do not 
contribute to the vv final state since neutrinos do not couple to Higgs bosons. The 
if -pole also gives a visible enhancement to the WW, ZZ and hh final states but not 
to W ± H T , ZA, HH, Hh, AA and H + H~ which do not become kinematically allowed 
until m x > ran/ 2. Likewise, the A-pole shows a resonance in the Zh final state but 
not in W ± H^, ZH, AH and Ah. 

For the pseudoscalar Higgs exchange in XX ~^ //> the expansion method gives 
a discontinuity very close to the A-pole, m x ~ = 200 GeV. For this channel, 

which is clearly dominant in the resonance region, the ^-coefficient is large and negative 
and the expansion method and iterarative procedure break down. This can be clearly 
seen in all the // windows (including vv) and also in the total value of J e xp(^/)- In 
the close vicinity of the poles, one can see the well-known large difference [0,[T^,|r7 



(note a log scale) between J exSuC t{xf) and J exp (xf) while away from the poles, both 
contributions are similar. 

We can also notice that J ex act(a ; /) an d J cxp (^/) show a significantly different en- 
hancement in the total values and in the up-type quark-pair final state in the region 
m x > m t = 175 GeV where the ti final state becomes kinematically allowed. In 
particular, J exact {xf) increases gradually starting from m x some 10 GeV below m t 
while Jexp{%f) exhibits a sharp increase only at m x = m t . This is because the exact 
treatment takes into account the thermal distribution of momenta of the annihilating 



neutralinos, while the expansion method neglects it [14]. The same effect is visible 
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near the thresholds of all the final states which only become open at some non-zero 
value of m x , especially of those involving one or two Higgs bosons. 

As we stated in Section |3], once kinematically allowed, the contributions from the 
W ± H^ final state can be comparable, or even larger (if sfermion masses are increased) 
than the // final state. In fact, the ZH, Ah and, at much larger m x also AH, final 
state contributions are not much smaller. We can also see that the channels ZZ and 
WW, as well as Zh, despite not being s-wave suppressed, are usually subdominant, 
except near m x ~ m A /2 because of the pseudoscalar Higgs exchange contributing to 
those channels, and when the higgsino admixture of the neutralino increases which is 
the case with increasing m x while keeping fi fixed. 

As a way of verifying our results, we have performed a numerical comparison of 
our exact results for the cross section (more precisely of the function w(s)) with the 
same quantity computed using DarkSusy. We have found an excellent agreement at 
the level of a few per cent for all the tree-level final states when we fixed the same 
values of input parameters. This gives us some confidence that the (highly complex) 
expressions that we have derived (and also those in DarkSusy), are correct. 

As mentioned above, in the current version of our numerical code we have used the 
iterative method (§) which is widely used in the literature. We applied it to both the 
exact formula for (avu0i), Eq. (|7|), and to the usual approximation fll3|) . The results 
are presented in Fig. [3] as a solid and dashed line, respectively. For comparison, we 
also display Xf using DarkSusy (dotted line) where a different numerical procedure is 
used. It is clear that there is a general agreement among the three different procedures 
and also that Xf can be only roughly approximated by the usually quoted value 0.05. 
Away from the resonances and thresholds the expansion method can be safely used 
in determining Xf while near such special points this is not the case. In particular, 
near the A-pole the expansion gives (av M01) < and the iterative procedure breaks 
down. In contrast, the exact numerical integration (|B|) provides reliable results both 
near and away from special points. However, we have carefully verified that this is the 
case only when the integration in Eq. @ is performed with high enough numerical 
precision. 

In Fig. |]a, we plot the ratio fi ex p/^exact versus m x for the same set of parameters as 
before. We denote by fl ex3lC t and f2 exp the neutralino relic density calculated with the 
exact (0) and approximate ([13]) formulae for the thermal average, respectively. The 
relic abundance in both cases is computed by solving Eq. (||) iteratively and using 
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Eq. From now on we impose all applicable experimental constraints, including 
the lightest Higgs mass constraint > 113.5 GeV |38|] and a lower bound 104 GeV 
for chargino masses. In the case of b — > 37 we have followed DarkSusy in using 
a somewhat old range 1.0 x 10~ 4 < B(b — ► S7) < 4.0 x 10~ 4 |39|] and have not 
yet included new SUSY contributions at large tan (3 that have been recently pointed 



out ||40|| . We feel that this is sufficient for the purpose of this analysis. The solid line 
represents the experimentally allowed region while the dotted line corresponds to the 
low m x range excluded by collider experiments. 

It is evident that the ratio f2 exp /f2 exact varies considerably with m x . Near reso- 
nances it is very large as discussed above. If we conservatively demand a 10% accuracy 
in the relic density computations then the approximate method fails to satisfy this 
criterion over a broad range 160 GeV < m x < 220 GeV, which is a significant fraction 
of the neutralino mass range. (In fact, various naturalness constraints, even if only 
indicative, give a rough constraint m x < 200 GeV |^|,|23|].) For larger values of tan/3 
the range of m x where 0.9 < ^ e xp/^exact < 11 is not satisfied becomes even larger. 

It is also clear that even further away from resonances the ratio f2 exp /f2 exact ex- 
hibits a rather complex behavior. For 220 GeV < m x < 260 GeV the new channels 
W ± H T , ZA and ZH° successively kick in, and become actually even somewhat more 
important than //. Since, as we mentioned above, the expansion method does not 



work properly near new channels thresholds |T4| , the ratio varies quite rapidly over the 
mentioned range of m x , although not as much as near resonances. Finally, far away 
from resonances and thresholds, the expansion method seems to work remarkably 
well. 

In Fig. |]b we show the relic abundance Q x h 2 as a function of m x . Again, the 
solid (dotted) line represents the range of m x allowed (excluded) by the experimental 
constraints. The two horizontal dashed lines delimit the cosmologically expected 
range ([]]) of Q x h 2 . Two features should be mentioned. Over the range of smaller m x , 
in particular near resonances, the relic abundance varies rapidly. On the other hand, 
at larger m x , where new channels involving the "heavy" Higgs bosons H ± , A and H 
open up and can become important, Q x h 2 exhibits a wide plateau. In both regions it 
is clearly important to compute Q x h 2 accurately if one wants to achieve an accuracy 
expected from future determinations of QcDwih 2 - In particular, it is clear from Fig. |4]b 
that at larger m x a relatively small difference of a few per cent in determining Q x h 2 
could exclude on cosmological grounds the range of m x > 280 GeV which in Fig. f|b is 
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(barely) allowed. This again illustrates the uncertainties involved in computing Q x h 2 
and in determining the cosmologically allowed regions of SUSY parameters. 

5 Conclusions 

While we are certainly still some time away from a "precision era" for measuring 
cosmological parameters, much effort has been focussed on improving theoretical cal- 
culations of the relic abundance of the neutralino. We have derived a full set of 
analytic expressions for the neutralino annihilation cross sections into all tree-level fi- 
nal states. In this paper, we have presented one example of such an expression for the 
neutralino annihilation into W + H~ which can be even dominant. We have also com- 
puted the first two terms in the expansion of the exact formulae in powers of x and we 
have compared the integral of the thermally averaged product of the cross section and 
velocity and the relic abundance computed both ways. We have confirmed the well- 
known inaccuracy of the approximate method near reasonances and threshold and 
showed that, far away from such cases the expansion method works well. However, 
due to the presence of several resonances and thresholds in the MSSM, the expansion 
method may lead to sizeable errors over a relatively large range of the neutralino mass 
of up to several tens of GeV. We have also demonstrated that the iterative way of 
computing the freeze-out point works rather well for both the exact and expansion 
methods of computing the thermal average but not in the latter case near resonances 
and thresholds. 
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Appendix 

Here we give expressions for the auxiliary functions used in the text. First, we define 



D(s,x,y 1 ,y 2 ) = x + 



yi + V2 _ s 

2 2 



(23) 




s 



(24) 
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If we define D = D(s, x, y 1 ,y 2 ),F = F(s, x, y u y 2 ), t±(s, x, y u y 2 ) = D±F and (%, y t ) 
= (T i ,y i )(s,x,y 1 ,y 2 ,z 1 , z 2 ) (i = 0, 1,2), then 



F{s, x, y h y 2 , z) = — In 



t+(s,x,yi,y 2 ) - z 



t_(s,x,y 1 ,y 2 ) - z 



(25) 



and 



1 



z\ - z 2 
1 

z\ - z 2 



r 2 = i + 



i 



z 1 + z 2 - 2 D 
1 



3>2 



zi + z 2 - 2 D 
1 



[F(s, x, 2/i, ?/2, «i) - ^(s, ac, 2/i, 2/2, 22)], 
[zi .F(s, x, 2/1, 2/2, *i) - ^2 F(s, x, 2/1, 2/2, z 2 )], 
\z\ T{s, x, 2/1, 2/2, zx) - z\ T{s, x, 2/1, 2/2, ^2)], 
[.F(s, a;, 2/1, 2/2, zi) + x, 2/1, 2/2, 2fe)], 
[2 (z x - -D) .F(s, x, 2/1, 2/2, Zi) - 2 (z 2 - D) F{s, x, 2/1, 2/2, 22)], 



Zi - z 2 
1 



1 + 



zi + z 2 - 2 D 



[zi fa - 2 D)F(s, x, 2/1, 2/2, 21) 



+ z 2 (z 2 -2D) F{s, x, 2/1, 2/2, z 2 )]. 



The expressions for Gjy% and G^jj, where i = 1, 2, 3 and j = 1, • • • , 5, in Eqs. (|17l - |19D 
are given by 



= s(m 2 +2m 2 v )-2m x + m 2 (m 2 H ±-m 2 v )-2m 2 v (m 2 H ±+m 2 



WJi 



G wh = (m 2 -m 2 H± )(m 2 -m 2 w )(m 2 +2m 2 w ), 



r Y(l) r T(l) 



WH i 



Gtfrf = 3 s m 2 w - 12 m 2 w m 2 x + 3 m 2 w m 2 H± - 3 m%r, 

G$§ = s 2 -s(2m 2 x + 2m 2 H± + 3m 2 w ) + 2mJ + (2m\ + m 2 H± )(m 2 H± + 3m 2 w ) -2m 
g wh = s(m 2 -2m 2 v )(m 2 -m 2 H ±)+Am 2 v m A 
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+ m 2 x (m 4 H± - 5m 2 w m 2 H± + 2m^) - 2m^ v m 2 H± , 
g wh = s m x^ m x _ m w) - m x ~ m^(m 2 H ± +3m 2 v ) + m 2 x m 2 v (2m 2 H ± + 3771^). 
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Figure 1: The total value of J(xf), Eq. (|^), and several partial contribution are 
shown in separate windows as a function of m x for tan (3 = 10, itlq = uiq = mu = 
m D = m L = mE = 500 GeV, m A = 400 GeV, A t = A b = 600 GeV and fi = -500 GeV. 
The solid lines represent the exact results, while the dotted ones correspond to the 
expansion (O. Notice that the final states W ± H Zf , ZH° and Ah, once kinematically 
allowed, give comparable contributions to the // channels. 
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tan/8 = 10, m = 500 GeV, m A = 400 GeV, A, = A b = 600 GeV, fi = -500 GeV 
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Figure 2: The same as in Fig. [I] but for mostly subdominant channels. Notice that 
in the lower two rows the horizontal axis has been shifted by 100 GeV. 
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Figure 3: The freeze-out point Xf as a function of m x . The solid and dashed lines 
corresponds to the iterative procedure (||) with {<jvm0i) computed exactly (|j) and in 
terms of the expansion flT3|) , respectively. For comparison, the dotted line has been 
obtained using DarkSusy. 
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tan/3 = 10, m = 500 GeV, m A = 400 GeV, A, = A b = 600 GeV, fj. = -500 GeV 
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Figure 4: The ratio ^exp/^exact (a) and the relic density Q x h 2 (b) for the same 
choice of parameters as in Fig. [I]. The solid (dotted) curves are allowed (excluded) by 
current experimental constraints. In window (a) the relic aboundance in both cases 
is computed by solving Eq. (|3|) iteratively and using Eq. (^). In window (b) we show 
Q x h 2 is computed using our numerical code. The band between the two horizontal 
dashed lines corresponds to the cosmologically favoured range 0.1 < fi x h 2 < 0.3. 
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